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We investigate the stability and geometrically non-linear dynamics of 
slender rods made of a linear isotropic poroelastic material. Dimensional 
reduction leads to the evolution equation for the shape of the poroelastica 
where, in addition to the usual terms for the bending of an elastic rod, we 
find a term that arises from fluid-solid interaction. Using the poroelastica 
equation as a starting point, we consider the load controlled and displace- 
ment controlled planar buckling of a slender rod, as well as the closely related 
instabilities of a rod subject to twisting moments and compression when em- 
bedded in an elastic medium. This work has applications to the active and 
passive mechanics of thin filaments and sheets made from gels, plant organs 
such as stems, roots and leaves, sponges, cartilage layers and bones. 

1 Introduction 

Poroelasticity is the continuum theory used to describe the behaviour of a 
biphasic material in which fluid flow is coupled to the elastic deformation of a 
solid skeleton (see Selvadurai (1996), Wang (2000) and references therein). 
The first applications of this theory were to geological problems such as 
consolidation of saturated soil under a uniform load (Biot 1941). Since then 
the theory has grown to cover many and varied applications, some of which 
are displayed in Table 1. 
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If a medium having interstitial fluid of viscosity v is forced to oscillate 
with a characteristic time r, the Stokes' length of the motion, L s = y/ur, will 
characterise the range of influence of the solid into the fluid. If Lg <C l p (the 
pore length scale) the fluid within the pores moves out of phase relative to 
the solid. On the other hand, if L s 3> l p , the fluid will only move relative to 
the solid when the volume fraction of solid matrix changes locally. This limit 
(L s S> lp) was first considered by Biot (1941) for an isotropic poroelastic 
material. Later work using averaging techniques led to equations of the same 
form as well an understanding of how the microstructure of the material 
influences the constitutive equations of the material (Auriault & Sanchez- 
Palencia 1977, Burridge & Keller 1981, Mei & Auriault 1989, Lydzba & 
Shao 2000). 

Whereas geological applications are concerned primarily with bulk be- 
haviour, many engineering, physical and biological applications have ex- 
treme geometries which allow for the application of asymptotic methods to 
reduce the dimension of the problem. Some examples include the active and 
passive mechanics of thin filaments and sheets made of gels, plant organs 
such as stems, roots and leaves, sponges, cartilage layers, bones etc. 

In this paper we use the constitutive behaviour of a linear isotropic 
poroelastic solid to investigate the stability and dynamics of slender rods 
made of this material. In £j2 we give a physically motivated derivation of 
the constitutive equations for a poroelastic material. In ^21 we use the bulk 
poroelastic constitutive equations to determine the equation for the time 
dependent bending of a slender poroelastic rod subject to an externally ap- 
plied compressive force P. Dimensional reduction leads to the equation for 
the poroelastica, where in addition to the usual terms due to the bending 
of an elastic rod we find a term due to the fluid resistance. It arises from 
the fluid-solid interaction and has a form similar to that in a Maxwell fluid 
(Bird, Armstrong & Hassager 1987). In 21 we solve the problem of load 
controlled buckling. Although the poroelastic nature of the material does 
not change the buckling threshold or the final stable shape, it governs the 
dynamics of the system as it evolves from the unstable to the stable state. 
Both the short and the long time limit are investigated investigated using 
asymptotic methods. We then use numerical methods to corroborate our 
asymptotic approaches and follow the non-linear evolution of the poroelas- 
tica. In we treat the problem of displacement controlled buckling and 
compare the results of the poroelastica with those of the classical elastica 
under similar loading conditions. In ^Hlwe consider the linear instability of 
a slender poroelastic filament embedded in an infinite elastic medium, when 
it is subjected to an axial twisting moment and an axial thrust. Finally, 
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Table 1: Chart illustrating various applications of poroelasticity. L s = \fur 
is the Stokes' length, where v is the kinematic viscosity of the interstitial 
liquid, and r is the time scale of the motion. 



in §7, we summarise our results and discuss possible applications to such 
problems as the mechanics of cartilaginous joints and rapid movements in 
plants. 

2 Governing equations for poroelastic media 

We begin with the equations for a homogeneous, elastic isotropic poroelastic 
material in the limit where the Stokes' length, L s = y/vr is much larger than 
the system size l m and further that l m is much larger than the pore size l p . 
We will also neglect inertial effects in the solid and liquid phases. In this 
limit the viscous resistance to fluid flow in the pores is balanced by the 
pressure gradient so that the momentum balance in the fluid yields 

^Vjv - Vp - V lpPp = 0, (1) 

where v is the fluid velocity with characteristic scale V, p is the fluid den- 
sity, V and Vi p denote gradients on the system scale and the pore scale 
respectively, p is the macroscopic pressure driving the flow, and p p is the 
microscopic pressure in the pore. When the pore scale and system size are 
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well separated (l p /l m *C 1) (2.1) gives the following scaling relations 

l m pvV pvV 

P w- » Pv (2) 

Thus the dominant contribution to the fluid stress in the medium arises from 
the pressure. The simplest stress-strain law for the composite medium then 
arises by considering the linear superposition of the dominant components 
of the fluid and solid stress tensor. Assuming that the elastic behaviour 
of the solid skeleton is well characterised by Hookean elasticity (i.e. the 
strains are small), we can write the following constitutive equation for the 
poroelastic medium (see Appendix A for a derivation using the method of 
multiple scales): 

a = 2 fie + AV • u I - apl. (3) 

Here a is the stress tensor, u is the displacement field, e = (Vu + Vu T )/2 
is the linearised strain, p and A are the effective Lame coefficients of the 
material (dependent on the material properties and the microstructure) , a 
is related to the fluid volume fraction, but includes a contribution from the 
pressure in the surrounding fluid (see Appendix A), and I is the identity 
tensor in 3 dimensions. These material parameters can be derived using 
microstructural information (see Appendix A). In the limit when inertia 
can be neglected, the equations of equilibrium are 

V • a = 0. (4) 

Mass conservation and continuity requires that the rate of dilatation of the 
solid is balanced by the differential motion between the solid and fluid in 
a poroelastic solid. This yields (see Appendix A for a derivation using the 
method of multiple scales) 

V • k • Vp = j3d t p + ad t V • u, (5) 

where the solid skeleton is composed of a material with bulk modulus (3~ l 
(7^ \ + 2p/3 since the lame coefficients A and p are for the composite ma- 
terial and take into account the microstructure, while is independent 
of the microstructure) and k is the fluid permeability tensor of the solid 
matrix. In words, © states that the flux of fluid into a material element 
is balanced by the change in solid volume due to the bulk compressibility 
of the matrix. For a rigid incompressible skeleton, /3 = 0. We will assume 
that the solid skeleton is isotropic so that k = fcl; however, many structured 
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Figure 1: Schematic diagram of (a) a bent rod, where 6{x) is the angle 
between the deformed and undeformed tangent vector, x, y and z are body- 
fixed coordinates in the reference frame of the rod; (b) the circular cross 
section. 



and biological materials are anisotropic and one may need to revisit this as- 
sumption. Equations (j3J), Q and © when subject to appropriate boundary 
conditions describe the evolution of displacements u and fluid pressure p in 
a poroelastic medium. The typical values of the parameters for a soft gel 
are a ~ 1, /x ~ A ~ 10 6 Pa, and k ~ 10 -12 m? /sec Pa. 

3 Equations of motion for a slender filament 

We consider a naturally straight slender circular rod of length L and radius 
R -C L with a tangent to the centre line that makes an angle 6 with the 
horizontal (see Figure At its ends an axial force P is applied suddenly 
at time t = 0. We will further assume that the lateral surfaces of the 
filament are free of tractions. This assumption could break down when the 
solid matrix is very dilute, so that interfacial forces become comparable to 
the internal forces in the filament, but we will not consider this case here. 
The slenderness of the filament implies that the axial stresses vary rapidly 
across the cross-section and much more slowly along it, so that we can use 
an averaging procedure to deduce low-dimensional equations that describe 
the motion of the filament. This long- wavelength approximation can be 
formalised using an asymptotic expansion in the aspect ratio of the filament 
R/L <C 1. Here, we will proceed directly by noting that since the rod is 
slender, bending it is easier than stretching or shearing it (Love 1944). At the 
level of scaling, geometry implies that the out-of-line (bending) displacement 
of the centre line scales as R while the axial displacement scales as R 2 /L. 
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At the surface of the filament no stress is applied. Since the filament 
is thin this implies that a yy ~ a zz ~ 0. For a displacement field u = 
(u x ,u y ,u z ) equation (J3J) yields 

a yy = ~ a P + ( 2 M + x )dy u y + K^xUx + <9 2 ii 2 ) = 0, (6) 
o" 22 = -op + (2/i + A)<9 2 u 2 + X(d x u x + = 0, (7) 

which can be solved for d y u y and d z u z to give 

a a ap-Xd x u x 

d y u y = d z u z = 2(/j + A) . (8) 

Equations © and (jSJ) give the axial stress 

. 3A/i + 2fi 2 

&xx — * ' P~\ OxU x . (9) 

A + fj, A + /i 

More specifically, when an infinitesimal axial element of the rod of length 
dx is bent so that locally the centreline curvature is d x 6, fibres that are 
parallel to the neutral axis (coincident with the centre line for a homogeneous 
circular cross-section) and at a perpendicular distance y from the neutral 
plane (defined by the neutral axis and the axis of bending) will be either 
extended or contracted by an amount yd x 9dx, so that the elastic strain 
d x u x = —yd x 9. This leads to an elastic stress that varies linearly across the 
cross-section; in addition there is a fluid pressure that is determined by ©• 
We insert (JHJ) into © and find the evolution equation for the fluid pressure 

k(d X xP + d yy p + d zz p) = ((3+ — — )d t p - -^y—yd xt 0. (10) 

A + [X A + jU 

To make the equations dimensionless we use the following definitions for 
the dimensionless primed variables: 

x = Lx', y = Ry , z = Rz', 6 = —6', 

Li 

(2fi 2 + 3Xfi)R 2 , vr(2^ 2 + 3\fi)R 4 



p = -ij; — -JiL — p> 



(n + A)L 2 xx ' 4(/i + A)L 2 

a^R 2 , [(fi + \)(3 + a 2 ]R 2 , ( 

P + x)/3 + a 2 ]L 2 P ' (ji + \)k { ' 

and the immediately drop the primes, referring exclusively to dimensionless 
variables from now on. We note that the axial stress, the compressive force 
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Figure 2: Schematic of the local torque balance in a bent rod under an 
externally applied compression P. We balance torques about the point O, 
x, y and z are the coordinates in the body-fixed frame with O as the origin. 
Let M(x) denote the total moment due to internal stresses generated in 
part by bending the elastic skeleton and in part by the fluid pressure field. 
Balancing torques gives: M(x + dl) — M(x) + dl x P = 0. In the limit 
dl 0, d x M + Psinf9 = 0. The total moment M = E p Id x 9 + a JpydA, 
where E p is the effective elastic modulus of the poroelastic skeleton and I 
is the moment of inertia of the cross section. 



and the pressure are scaled to reflect the dominance of bending deformations 
over all other modes, and the time is scaled to reflect the dominance of radial 
diffusion over axial diffusion. 

Then the stress in the filament, given by equation @, can be written in 
dimensionless form as 

cr X x = -y d x 6 - -p (12) 

Here, the first term reflects the purely elastic contribution well known from 
the theory of beams (Love; 1944), while the second term is proportional to 
the fluid pressure in the pores. The dimensionless parameter 5 



0(1) for most materials denotes the ratio of the fluid and solid stress. 

In the long wavelength approximation, it is preferable to use the stress re- 
sultant F = f o xx dA and the torque resultant M = — f ya xx dA = f y 2 d x 9dA+ 
j f ypdA = M e + Mf, as the variables of interest. Here M e is the elastic 
torque and Mf is the fluid torque that arises due the transient effects of a 
pressure gradient across the filament. Then, local force and torque balances, 
which can be derived from ©, (J3J) and ©, or equivalently directly (Figure 
El, yield the dimensionless equations 



d x F = 

p 



d x M + ^P sin 9 = 0. (13) 
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The first of these equations can be integrated immediately to yield F = P 
with P a constant determined by the boundary conditions. The second 
equation combines the effects of the elastic and fluid stresses that arise due 
to the fluid pressure, and requires the solution of the continuity equation 
(|1U|) . For a rod with a circular cross-section, there is rotational symmetry in 
the problem. Choosing the the axis of bending to coincide with the z axis, 
we rewrite in polar coordinates (r, <fi) using dimensionless variables as 

d t p - -d r (rd r p) - -zdssp = r sin (/)d xt 9. (14) 

We see that the pressure in the fluid arises from the extensional and com- 
pressional stresses in the filament due to bending. The boundary conditions 
for the pressure can be deduced using the following considerations: (a) the 
centre line of the rod does not suffer any deformation, and is symmetrically 
disposed, and (b) the pressure at the surface is determined by the perme- 
ability of the surface layers and the flux through it. Then 

p = at r = 0, 
Bip + d r p = at r = l. (15) 

where Bi = and 77 characterizes the flux through the surface for a given 
pressure drop (the ambient external pressure is assumed to be zero). The 
second boundary condition in IJ15JI on the pressure states that the flux of 
fluid through the surface is proportional to the pressure drop across the 
surface. Bi = oo corresponds to a freely draining rod, where there is no 
pressure jump across the surface, and Bi = corresponds to a jacketed rod, 
which allows no flux through the surface. For a sponge, Bi> 1, while for a 
plant (root, stem or leaf) Bi < 1 since it is designed to retain water. 
Expanding p in terms of the homogeneous solutions of (|14|) . we write 

oo oo 

P = ^2 ^2 [A mn sin mcj) + B mn cos m<p] J m (r V / A^)e- Am "*, (16) 

m=0 n=l 

where A mn and B mn are constants, J m is the Bessel function of order m, and 
Xmn is determined by the boundary conditions. Inspection of the inhomo- 
geneous term on the RHS of equation (fTl*|) yields m = 1 so that A\ n = A n , 
B mn = 0, and Ai n = A n . Since the boundary condition is a linear combina- 
tion of p and d r p we are guaranteed to have a complete basis. We therefore 
look for a solution to the inhomogeneous equation (|14|) of the form 

oo 

p = ^2A n {t)shx4>Ji{r^/Xn), (17) 

n=l 
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where A n is determined by substituting 1)17(1 into lfTo)l which yields 

d r Ji{r^JJ^) +Bi Ji(r^/\^) = at r = 1. (18) 
Inserting (|17|) into (fTI)) yields 



^2(d t A n + A^J^r^) = rd xt 9. (19) 



n=l 



Multiplying (|19|1 by r J\{r^J\ n i) and integrating across the cross-section 
gives 

d t A n + X n A n = Xn d xt 9, (20) 

where 

Jor[Ji(rv^)] 2 *' 

Solving equation (|2U() yields 

A n = Xn t e-^-^d^Odt', (22) 

JO 

so that (|17() may be rewritten as: 

oo „ t 
p = Y i XnSixi<f>Ji(rVK) e-^-^d^Odt'. (23) 

n=l ^ 

Then (|12() allows us to write the total axial stress a xx at a cross-section as 
cr ra = -r sin 03,0 - - J2xn sin M(rVA^) / e~ A ™(*-* )d xt ,9dt'. (24) 

n=l ^° 

The dimensionless torque resultant is given by 

r sin <jxj xx dA = -d x 9 + ^ r Y J 7n e~ Xn ^ d xV 9dt' , (25) 
4 4 „=i ^ 

where 7„ = Xn J r 2 Ji{r\fX^)dr . Substituting the result into the equation 
for torque balance (|13l) yields the dimensionless equation for the poroelastica 
(see Figure HJ) 



oo ,. t 

d xx 9 + Psm9 + 8Y jln e- x ^-^d xxtl 9di! = 0, 



(26) 



Fluid Resistance 



HMIN 



Load <- 




Load 



Elastic Resistance 

Figure 3: Mechanical analogue of the bending resistance of a poroelastic rod. 
For rapid displacements, the dashpot will not move and the fluid resistance 
due to the instantaneous pressure yields a response similar to a stiff (fluid) 
spring in parallel with an elastic spring. Eventually the dashpot will move 
to relieve the stress in the spring and the fluid resistance gradually decays, 
leading to a purely elastic steady state. 



The first two terms correspond to the usual terms in the classical elastica 
(Love 1944) for the bending of a rod with a circular cross section, while the 
final term is due to the instantaneous fluid pressure not being equilibrated 
across the cross section. The influence of the fluid is to create a material 
with "memory" , so that the current state of the filament is determined by 
its entire history. The kernel in the memory function for the fluid resistance 
is e~ x "( t ~ t ) so that the fluid resistance is analogous to that of a Maxwell 
fluid (Bird, Armstrong and Hassager 1987) with relaxation times 1/A ra which 
measure the rate of decay of the n th transverse mode in response to the rate 
of change of the curvature of the filament d x t9. A mechanical analogue of 
the resistance of a poroelastic filament is presented in Figure Eland shows 
the connection to simple viscoelastic models. 

The dynamics of a poroelastic rod are then determined by the solution 
of the equation for local torque balance (|26|). subject to the boundary con- 
dition on the fluid pressure at the surface which determines the decay 
constants A n , and additional boundary conditions on the ends of the rod, 
which we now consider in some specific cases. 
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4 Planar load controlled buckling 



When an initially straight rod that is simply supported at either end is 
subject to a constant compressive force P applied suddenly at t = 0, the 
boundary conditions at the ends are 

0*0(0, t) = o, 

d x 6(l,t) = 0, (27) 

and the initial condition is 

9(x,0) = 0. (28) 

The complete time evolution of the rod is then given by the solution 
of the integrodifferential equation ()26|) subject to the boundary conditions 
(|27|). the initial conditions (|28(1 and the condition (|18|) which determines the 
rate constants A„. 



4.1 Short time behaviour, t<l 

Expanding the solution about the initially straight state 9 = 0, we write 

9 = e9 x (t) (29) 

where e« 1. Substituting the expression (|29|) into equation (|26|) and lin- 
earizing yields 

oo „ t 

d xx 9 x + P9 1 + 5Y jln e- x ^')d xxt ,9 x di! = 0, (30) 

n=l ^0 

subject to the boundary conditions 

d x 9 1 (0,t) = d x 9 1 (l,t)=0. (31) 

To solve ()30II31|) we use separation of variables writing 9±(x, t) = g(x)f(t) 
and substituting the result into (|30|) to obtain two equations for g(x) and 
f(t). The function g(x) is determined by the solution of the eigenvalue 
problem : 

(l + 09 xx g + Pg = 0, d x g(0) = 0, d x g(l)=0. (32) 

Here the separation constant £ = [P — ir 2 )/ir 2 is the relative difference 
between the applied load P and the dimensionless buckling load, P c = it 2 
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for a purely elastic rod that is simply supported at its ends. The function 
f(t) satisfies 

= X> fe-^-^d.mdt', (33) 

Using Laplace transforms (C(f(t)) = / °° e~ st f(t) dt) we solve (|3*3*)) and find 
that 

oo 

/(t) = -£*•*/(<)) (34) 
seS n=l An + s 

where /(0) is determined by the initial condition and the set S is composed 
of elements which satisfy 

& oo 

The growth rate at the onset of the dynamic buckling instability is therefore 
given by the largest s that satisfies (|35|l . In Figure 0J we plot the growth 
rate s as a function of the rescaled external load | = P ^ c , with P c = ir 2 , 
obtained by solving (|32|) and ()35|) . When s < we do not have an instability, 
corresponding to the case when P < P c . In the poroelastic regime, when 
P c < P < -P c (l + |), fluid flows across the filament in response to the stress 
gradient in the transverse direction, and the phenomenon is qualitatively 
different from the buckling of a purely elastic rod. Since the time it takes 
fluid to flow across the filament is longer than it takes a bending wave to 
propagate the length of the filament, poroelastic buckling is sometimes called 
creep buckling (Biot 1964). 

We now turn to the dependence of the buckling transition on the sur- 
face permeability parameter Bi. Substituting a pressure field of the form 
p(x,r,(f),t) = h(r)d x Q(x) sin <pe st and 9 = e st Q(x) into Q14JI yields 

sh — -d r (rd r h) + 4j = sr, (36) 

subject to boundary conditions (J15|) which are now given by 

h = at r = 0, 
Bih + d r h = at r = 1, (37) 

Thus (|3*5|l yields lim^o h = and lim^oo h = —r corresponding to the case 
of infinitely slow and infinitely fast growth rates respectively. Consequently, 
for infinitely slow buckling the fluid supports no load. In the case s — > oo 
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Figure 4: Growth rate s of the deflection as a function of the dimensionless 
external load P; P c = tt 2 is the critical compression above which a simply 
supported purely elastic rod buckles. Here the surface permeability param- 
eter Bi= 0.1. When P ^T C = \ the growth rate becomes infinite and one 
must consider inertial effects, which are neglected here. For comparison, we 
show the numerical results obtained by solving Q26JI with initial conditions 
9(0) = 0.001 cos 7tx, 9(dt) = 0.001 e sdt cos ttx, where s is the theoretically 
calculated exponent, with dt = 0.001, dx = 0.01. 
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Figure 5: The radial variation of the fluid pressure at the onset of buckling, 
h(r), (a) for growth rate s = 1 and various values of the surface permeability 
parameter Bi (larger Bi corresponds to a more permeable surface) and (b) 
for Bi = oo and various s. 



near r = 1 a boundary layer emerges where the internal solution (h = r) 
is matched to the boundary condition (|37|) at r = 1. Balancing the first 
two terms of (|36(l the length scale of the boundary layer Z&j ~ or in 

dimensional terms Li ~ \/ -n ^tt^t si • 

To complement these asymptotic results, we solve ([3613 7j) numerically 
and plot the radial variation in the pressure. In Figure we show h(r) 
for s = 1 and various Bi and in Figure we show h(r) for Bi = oo and 
various s. As expected, we see that as the surface permeability increases 
(i.e. Bi increases) for a given growth rate of the instability (corresponding 
to a given load) the pressure variations across the filament decrease. On the 
other hand, as the growth rate increases, a boundary layer appears in the 
vicinity of the free surface of the rod to accommodate the slow permeation 
of fluid in response to the stress gradients. 

Having considered the onset of poroelastic buckling we now turn to the 
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transition from poroelastic to inertial dynamics which occurs for very large 
compressive loads when the fluid cannot move rapidly enough to keep up 
with the elastic deformations. For large s the condition that determines the 
growth rate l|3*5|) reads 

»!-$>«(!- V*)- (38) 

n 

We have computed En 7™ = 1/4 for integrals of Bessel functions. Using the 
definition of the separation constant £ = P/P c — 1 we solve equation (|38|) 
for the growth rate: 

„ _ En 7nA„ 



P-Pc 1 



> (39) 

8P C 4 

showing that it indeed diverges when F ^T C — * \ consistent with Figure 0] 
4.2 Long time dynamics t 3> 1 

In the long time limit t > 1, i.e. when the fluid has enough time to diffuse 
across and along the filament, the shape of the filament approaches that of 
the ideal elastica. To capture the dynamics of this process, we linearise (|26|) 
about the steady state solution by letting 

9 = 9 (x)+ee 1 (x,t) (40) 

with e <C 1. Substituting the expansion (|4*U|) into (|26|). at leading order we 
get 

d xx 9 + P sm6 = 0. (41) 

At 0(e), we get 

d xx e x + P cos(0 o )0i + <^7n /V A " cW Mi' = 0. (42) 

n ^ 

To simplify the equations further we consider the convolution integral in 
(|42jl for typical values of Bi = 0.1, corresponding to the case for soft gels 
and biological materials. Then |JT8|) yields {A n } = {3.67, 28.6, 73.1, 137, 
221, 325, ...} and Y^=2^nhi = 0.0152. Given the large separation between 
the decay constants, we see that the dominant contribution in the integral 
arises from Ai leading to an approximation of (|42j) that reads 



^ 1 + Pcos(, ), 1 + , 7l /e-M t -O^, M ^ . 



(43) 
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P-Pc 
6Pr 



Figure 6: As the solution approaches the equilibrium shape the difference 
between the current and equilibrium shape decays exponentially with a rate 
if), which is plotted against the applied compression. 5 = 1, Bi = 0.1, 
dx = 0.01, dt = 0.01. The numerical computation is begun with 9 = 0.9 ho- 



using separation of variables, 9\{x,t) = g(x)f(t), we find that g(x) is given 
by the solution of the eigenvalue problem 

(1 - Odxxg + P cos(0 o ) 9 = 0, d x g(0) = 0, d x g(l)=0, (44) 

while the temporal part f(t) satisfies 

~/£ = <57i f e- x ^-^d t ,fdt'. (45) 

Since we are interested in the asymptotic behaviour for t > 1 we multiply 
both sides of equation (|45|) by e Al * and differentiate with respect to time to 
find 

dtf = f^-J = -i>f (46) 

We observe that the poroelastic solution approaches the elastic steady shape 
exponentially fast at late times. In Figure |SJ we plot the exponent if) = 
(^7i + £ ) versus (P — P c )/ 5P C and see that the larger the value of P the 
faster the solution approaches the final shape. 

4.3 Intermediate time dynamics 

For intermediate times we have to solve for the shape of the poroelastica 
numerically. Our arguments in the previous section allow us to neglect the 
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contributions from the higher modes so that a good approximation to (|26|) 
is given by 

r-t 

8 XX 9 + Psin 9 + 7l / e- x ^-^d xxt ,ed£ = 0. (47) 
Jo 

For ease of solution, we convert the integrodifferential equation to a partial 
differential equation by multiplying Q47[) by e Al * and differentiating with 
respect to time, so that 

(1 + ji)d xxt 6 + Xxd xx 9 + P cos 9 d t 6 + Ai P sin = 0. (48) 

We solve equation (f!5|) subject to the boundary conditions (fTTj) using a 
Crank-Nicolson finite difference scheme in space and we extrapolate the 
non-linearity using the previous two time steps. This gives us a scheme 
with second order accuracy in time. For a time step dx = 0.01 and a space 
step dt = 0.001 the difference between the numerical and analytical initial 
growth rate is 0.2% (see Figure^}. In Figure [7|we show the variation of the 
angle 9(0, t) determined using the numerical simulation for the case when 
the dimensionless buckling load is slightly larger than the threshold for the 
poroelastic buckling, with (P — P c )/P c ~ 0.17. For comparison, we also 
show the asymptotic solutions for short and long times determined in the 
previous sections, and find that they agree well with the numerical solution. 
To determine the shape of the filament we use the kinematic relations 

d x X = cos 9, d x Y = sin 9, (49) 

where X(x,t) and Y(x,t) are the position of the centreline. Figure |H] shows 
the shape of the filament as it evolves from the initially unstable straight 
shape to the final elastic equilibrium via a transient overdamped route. In 
sharp contrast, a purely elastic rod subject to the same initial and boundary 
conditions would vibrate about the final state forever (in the absence of any 
damping) . 

5 Displacement controlled planar buckling 

In many problems involving instabilities, there is a qualitative difference 
between load controlled and displacement experiments. To understand the 
difference we consider the problem of displacement controlled buckling of 
a poroelastic filament and compare the results with those of the previous 
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Figure 7: 0(0, t) for P = 11, 6 = 1, Bi=0.1, dx = 0.01, dt = 0.001 and 
9(x,0) = 9(x,dt) = 0.01 cos irx. The short time asymptotic is for a growth 
rate s(P) found from equation ()35j) . The long time asymptotic is of the form 
A e -i>t _|_ where tp is the rate of decay to the equilibrium angle 9q, and A 
is a fitting parameter. 




o.o 0.4 0.8 

X 



Figure 8: Shape of the buckling filament X(t), Y(t) for P = 11, in the 
laboratory frame as a function of time for P = 11, 5 = 1, Bi=0.1, dx = 0.01, 
dt = 0.001 and 6(t = 0) = 0(i = dt) = 0.01 cos ttx. 
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section. Since the centre line of the filament is assumed to be inextensible, 
the change in the end-to-end distance is given by 



A(t) = 1 - I" 
Jo 



cos 6(x,t)dx. (50) 



We choose the functional form 



A(t) = =^[l + tanhot], (51) 

to allow us to ramp up the displacement to a maximum amplitude A max 
at a characteristic rate ^A max a. The shape of the poroelastica is now de- 
termined by the solution of l|48|l. (|5U|l and (|51|): the unknown load P(t) is 
now determined at every time step by using an iteration method to enforce 
(|5U|) . For an initial guess to start this procedure, we note that after the 
onset of buckling when P > P c , for small amplitudes 6 = e cos7rx (e « 1) 
is a solution of (|4"%|) and (|2Tj) . Substituting into lf5Ti|) gives 



f 1 e 2 
A = 1 — / cos(e cos7rx) dx — . 
Jo 4 



(52) 



Therefore, we choose 9(x,to) ~ 2^J A(to) cos irx. In Figure Eh, we show 
the evolution of the load P(t) for various values of a. P is roughly constant 
for very short and very long times, but changes as A varies quickly for 
intermediate times. We can understand the initial plateau by considering 
the case when (e at <C 1), so that (f5lj) yields 

A = %^[l-^S]-A_e- (53) 

In light of the geometrical constraint ()52|) valid for small displacements, this 

yields 

6 sn 2^A max e at cos ttx. (54) 

Comparing this with the short time behaviour of a poroelastic filament con- 
sidered in § 44. II we see that exponential growth of small angles corresponds 
to a constant compressive force seen in Figure 03 A similar argument holds 
for late times, when the system relaxes to its purely elastic equilibrium. 
For intermediate times, the load can be larger than that for the case of a 
purely elastic filament. The difference in the loads is due to the fluid resis- 
tance in the porelastica. A way of visualising this is shown in Figure ITU1 For 
slowly applied displacement fields P max is almost the same for the elastic and 
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Figure 9: P fjT c , where P is the load and P c is the critical load required 
for buckling, corresponding to a change in end to end displacement A(t) ~ 
0.1[l + tanh at], for various a. For some later times the more quickly applied 
displacement corresponding to larger a requires a lower compressive force. 
The graphs correspond to the following parameter values: dx = 0.01, dt = 
0.002, Bi=0.1, 6 = 1, 0(-3) = 0(-3 + dt) = 2 v / A(-3) cos irx. 
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poroelastic cases; however, for rapidly applied displacements, corresponding 
to large values of a, the compressive force in the poroelastic case is larger 
due to fluid resistance arising from the pressure gradients across the bending 
filament. 
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Figure 10: Maximum compressive load , P max , for displacement controlled 
buckling, where the displacement field is given by A(t) ~ 0.1[1 + tanh at]. 
Bi = 0.1, 5 = 1, dt = 0.01, dx = 0.01. 



6 Filament embedded in an external medium and 
subject to axial torque and axial thrust 

We finally turn to the case of a rod embedded in an external medium subject 
to an axial moment, K, and a compressive force P (see Figure ITT|). The 
presence of the twist causes the instability to become non-planar, and the 
filament adopts a helical conformation; the presence of an external medium 
typically causes the instability to manifest itself with a higher wave number 
than otherwise. Letting the displacements of the centre line in the y and z 
directions be Y(x,t), Z(x,t) respectively, we scale the kinematic variables 
accordingly to define the dimensionless displacements 

Y = RY', Z = RZ'. (55) 

The dimensionless axial moment is defined as K = -^-^^^—K'. If the 
transverse displacements of the filament are small, the resistance of the 
external medium can be well approximated using the response of a linear 
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F external 




Figure 11: Schematic diagram of a rod buckling under an applied twist and 
compression in an external medium. Y(x,t) and Z(x,t) are the displace- 
ments of the centre line in the y and z directions respectively. The dotted 
line denotes the axis of symmetry. 



Hookean solid. In light of the analogy between linear elasticity and Stokes 
flow, we can use the results of classical slender body theory in hydrodynamics 
(Batchelor 1970, Cox 1970) and write the vector of dimensionless external 
forces on the filament as F external = 

(0, -ttE Y/4, -ttE Z/4), where the 

dimensionless parameter E is given by 

= 16 Mm L 2 (/, + A) 

ln(^)^ 2 (2/i 2 + 3/xA)' 1 ! 

where \i m is the Lame coefficient of the surrounding medium. This approx- 
imation is valid when R/L <C 1, a condition consistent with the geometry 
of a thin filament. We will further assume that the filament is free to rotate 
in the medium, i.e. there is no torque resisting this mode of motion, which 
varies in any case as R 2 and is thus negligible in most situations. 

To derive the evolution equation for the shape of the filament we use the 
constitutive equation (J3J) to write down equations for the balance of forces in 
the y— and z— directions as for the planar filament. After dropping primes 
this leads to 

oo „ t 

d xxxx Y + Kd xxx Z + 5 V 7 n / e- x " {t -^d xxxxt Ydt' + Pd xx Y + £ Y = 0, (57) 
n=i J ° 

CO , t 

d xxxx Z - K8 XXX Y + iV 7 n / e~ A ™ W> d xxxxt Zdt' + P8 XX Z + E Z = 0, (58) 
n=l Jo 

where, since equation (fTH) is linear, we have superposed the two solutions 
for bending in the y— and z— directions. Taking £ = Y + i Z equations Q57|) 
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and (|58|) may be written as a single equation for the complex variable £ 



= d xxxx (-iKd xxx (+Pd xx (+E(+5J2'yn e-^-^d^^Cdt'. (59) 

n=l J 

For a simply supported filament, the four boundary conditions are 

C(0) = C(l) = d xx ((0) - i K8 X ((0) = d xx C(l) - % Kd x C(l) = 0. (60) 

We can treat equations (|59j) - H6U[) in exactly the same fashion as the planar 
problem and use separation of variables C( x ,t) = 9{ x )f{t) to get 

(1 + i)d xxxx g - i Kd xxx g + Pd xx g + Eg = 0, 
5 (0) = g(l) = d xx g{0) - i Kd x g{0) = d xx g(l) - i Kd x g(l) = 0, (61) 

an eigenvalue problem for the separation constant £ and g{x). The temporal 
part of the solution f(t) satisfies 

e oo 

|_V-^ = 0. (62) 

which is the same as equation (|45[) for the temporal part of the solution for 
planar buckling. Thus, once the separation constant £ is found, equation 
(|35|) yields the growth rate s(P, K, E, 5, Bi) as a function of the loading 
parameters and the material constants (see Figure 

As an example of how the influence of an external medium can lead to 
higher modes becoming unstable at lower compressions than the fundamen- 
tal mode we consider equation (|61|) in the case K = 

(I • • • AX 0. (63) 

which is an eigenvalue problem for £ and £ for a given P. At the ends 
[x = and x = 1) the displacements and bending moments vanish so that 
the boundary conditions associated with (|63j) are: £(0) = d xx C(0) = C(l) = 
d xx C(^) = 0. The only nonzero solutions to ()63|) occur when £ = sin q n x, 
where q n = rnr, n = 1,2,... The critical compression P c (n) where the n th 
mode becomes unstable is found to be (Landau & Lifshitz 1970) 

P c (n) =7rV + -=-=. (64) 

We see that the critical buckling load for a given mode number n increases 
as the stiffness of the environment E increases. Furthermore, for large E, 
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Figure 12: For the case of an applied compression P and axial twisting 
moment K, we show the three short time regimes (stable, poroelastic, and 
inertial) as in Figure 0] E = 5 = 1. 



the chosen mode shape does not correspond to the fundamental mode n = 1, 
since dP c /dn = yields n = ^E 1 ^ for an infinite rod. Physically this occurs 
because short wavelength modes do not deform the stiff elastic environment 
as much, while the penalty associated with a higher curvature is not too 
much of a price to pay. 

For the case when K ^ 0, we cannot solve the eigenvalue problem analyt- 
ically, and present the results using a phase diagram shown in Figure IT21 We 
find three distinct regimes: for s < the system is stable; for < s < oo we 
have the poroelastic regime where the system buckles on the same time scale 
as the fluid pressure diffuses; finally, we have the elastic regime where the 
system buckles so fast that the fluid does not move and all the deformation 
occurs in the solid skeleton. 



7 Discussion 

The usefulness of poroelastic theory is limited to a range of time scales. The 

poroelastic time scale associated with decay of pressure fields is t p ~ a ^ , 

recalling that a is the fluid volume fraction, R is the smallest macroscopic 

length scale of the system (R 3> l p ), is the effective Lame coefficient of 

i 2 

the composite material, and k ~ is the matrix permeability. If the time 
scale of the forcing, r <C t p the fluid will not move relative to the solid 
and Hookean elasticity and the effects of inertia are sufficient to describe 
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the system adequately. If r » t p , the fluid pressure will equilibrate with 
the surroundings and once again classical elasticity suffices to describe the 
system, albeit with different Lame coefficients. However, if r ~ t p the 
dynamics will be governed by poroelasticity. 

Biological systems are composed mainly of fluid, so poroelasticity will 
be applicable at some time scale (see Table 2 for estimates). Furthermore, 
they are characterized by extreme geometries (e.g. beams, plates and shells), 
which led us to consider in detail the dynamics of slender poroelastic objects, 
and particularly the buckling of a planar filament. Biological materials are 
usually anisotropic and we expect the permeability and elasticity tensors to 
reflect this feature. Taking ki to be the permeability in the axial direction, we 
can neglect axial diffusion if <C 1. The opposite limit, where S> 1, 
has been studied by Cederbaum et al. (2000). The dynamical behavior of 
these objects is separated into two different regimes, one governed by fast 
inertial effects, and the other by the slow dynamics of fluid flow. These 
regimes are of course well known in bulk materials (Biot 1957), but here 
they appear in a slightly different guise due to the effect of the slender 
geometry of the system. The onset of planar poroelastic load controlled 
buckling was first considered by Biot (1964). In this paper we broaden and 
deepen the understanding of this phenomenon. An important outcome of 
our studies is the poroelastica equation, which is a simple integro-differential 
equation with one time constant that describes the dynamics of a poroelastic 
filament under a compressive load. The bending resistance of the filament 
is analogous to a (fictional) Maxwell material, where the time constant is 
the rate at which the pressure field decays (determined by the material 
parameters and the geometry). We then used this equation to study not 
only the onset of buckling, but also the entire dynamics up until saturation 
for both load controlled and displacement controlled buckling. 

A series of three-point-bending experiments (Scherer 1992; Scherer 1996) 
have shown that the mechanical response of a silica gel rod immersed in 
acetone or ethanol can be described using poroelastic theory. The theory 
developed by Scherer (1992) applies only to situations where the displace- 
ment is applied much faster than the poroelastic time scale and is a special 
case of the more general theory presented in this paper. The lack of exper- 
iments on slender poroelastic filaments being deformed on the poroelastic 
time scale prevents us from testing our predictions quantitatively. Since our 
results are relevant to swollen polymer networks, gel actuators and sensors, 
the mechanics of cartilaginous joints, and the physics of rapid movements 
in plants, an important next step is the quantitative experimental study of 
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Application 


fi (Pa) 


a 


k (m 2 /Pa sec) 


R(m) 


t p (sec) 


Actin cytoskeleton 


100 


0.8 


io- 12 


10~ e 


io-' 2 


Bones 


10 iu 


0.05 


io- 14 - 10" ib 


10~ 2 


0.1 - 10" 3 


Cartilage 


10 e 


0.8 


(1 - 6)10 _ie 


10~ 3 


10 a 


Plant stem/root 


10 8 


0.8 


io- n 


10" 2 


10~ 2 


Venus' fly trap leaf 


10 e 


0.8 


io- 1 ^ 


io- 3 


0.1 



Table 2: Applications of poroelasticity in biology. 



slender poroelastic structures. 
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*Appendix A : Derivation of the poroelasticity equations 
The derivations of the equations of poroelasticity have been many and 
varied. Partly this has been because several qualitatively different parameter 
regimes containing distinct leading order force balances exist. We focus here 
on the equations which govern the second row of Table 1, namely where the 
Stokes' length is much larger than the pore size, i.e. L s 3> l p . The methods 
used to derive equations for this region of parameter space can be classified 
into three categories: physical arguments and superposition (Biot 1941, Biot 
& Willis 1957), mixture theory (Barry & Holmes 2001), and micro structural 
derivations (Auriault & Sanchez-Palencia 1977, Burridge & Keller 1981, Mei 
&; Auriault 1989). First, we show in detail a version of the micro structural 
derivations, which uses ideas from both Burridge & Keller (1981) and Mei 
& Auriault (1989). 

The equations that govern the behavior in the incompressible interstitial 
fluid at low Re are 

a f = - P I + 2e^e(v) , (A-65) 
V-<r/ = 0, (A-66) 
V • v = 0, (A-67) 

where <7j is the stress tensor in the fluid, e = l p /l m (see figure IT3|) . e(..) = 
tj[V(..) + V(..) T ] is the strain operator, and v is the fluid velocity. In the 




Figure 13: Typical porous medium illustrating the separation of length 
scales, where l p is the pore scale and l m is the system scale. 

solid the analogous equations are 

a s = A : e(u), (A-68) 
V • a s = 0, (A-69) 

where u is the displacement field, A is the tensor of elastic moduli, and a s 
is the stress tensor in the solid. At the solid-fluid interface continuity of 
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displacements and tractions yields 



v - d t u -- 
n — a j ■ n 



0, 

: 0. 



(A-70) 
(A-71) 



Here n is the unit normal vector to the surface separating the two phases. 

Looking for a perturbation solution in powers of the small parameter e, 
we use an asymptotic expansion of the variables 



°7 = 
cr s = 

P 

U : 
V 



<rj + e<r} + 

■ <T° S + £<r\ + 
= p° + ep 1 + 
= u° + eu 1 + 
= v° + ev 1 + 



with a multiple-scale expansion for the gradient 

V = X7 X , + eV, 



(A-72) 



(A-73) 



where x denotes the macroscopic scale, x' = ex denotes the pore scale, V 
denotes the gradient relative to the macroscopic scale and V x / denotes the 
gradient relative to the pore scale. Since we assume that the flow is driven 
on the macroscopic scale, the leading order deformation is a function only 
of x. Then equations ()A-65|) and (|A-68|) yield the following expressions for 
the fluid and solid stress tensors: 



A : [e(u°) 



V)], 



A : [efu 1 ) + e a ./(u 2 )], 

CT° f = - P °I, 

Vl + ^v(v°), 



(A-74) 
(A-75) 
(A-76) 
(A-77) 



where e and e x / denote the strain relative to the system scale and pore scale 
coordinates respectively. The stress balance in the fluid (|A-69|) yields: 



/,V 2 v° 



V x ,p l - Vp° 



0. 

= 0. 



(A-78) 
(A-79) 



Thus the leading order pressure gradient p°(x) is only a function of the 
system scale coordinate. The stress balance in the solid yields 



V, 



'i + v 



0. 

: 0. 



(A-80) 
(A-81) 
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we define a to be the total stress tensor: 

a = a s in V s , 

a = <Tf in Vf, (A-82) 

where V s and Vf are the solid and fluid parts of a volume element. Stress 
balance in the fluid and solid imply 

• a 1 + V • cr° = 0. (A-83) 

Averaging (|A-83|) over the pore scale. 

^ J v • a°dv + ^ y • o-W = 1 y v • + 1 y n • o^dS = 0, 

(A-84) 

where V = Vy + V^. In the limit F — ► oo, ^ J n • ^dS —> since the surface 
to volume ratio tends to zero. Consequently, 

i / V ■ er°dV = V- < cr° >= 0, (A-85) 



V _ 

< cr° > = < A : [e(u°) + e^u 1 )] > -<f> f p°I, (A-86) 

where 0/ is the fluid volume fraction and < > denotes averages over the pore 
scale. In order to write the averaged equations in terms of u° and p°, we 
must eliminate u 1 . This is achieved by using the stress balance in the solid 
so that 

V*' • a° s = V x , ■ {A : [e(u°) + e^u 1 )]} = 0. (A-87) 
The boundary condition (|A-71|) at the fluid solid surface yields 

A : [e(u°) + e x /(u 1 )] • n = -p°n. (A-88) 

Since this is a linear system of equations, u 1 is a linear combination of p° 
and e(u ): 

u 1 = B : e(u°) - Cp°, (A-89) 

where the third rank tensor B and vector C vary on the pore and system 
scales, and can only be found explicitly by solving the micro structural 
problem (T07l) - (fA^88l) . The averaged stress tensor becomes 

< cr° >=< A + A : e x /(B) >: e(u )- < A : e x /(C) > p° - <f> f p°I, (A-90) 

where in index notation e x /(B) = {d x i B n ki-\-d x i B m },i) If we assume that the 
material is isotropic on the macroscopic scale we can further reduce (|A-90f) : 

< a >= 2^e(u°) + AV • u°I + (-0/ + 7 )p°I, (A-91) 
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where 7I =< A : e(C) > is an isotropic pressure in the solid due to the 
fluid pressure exerted at the interface. Substituting (|A-91|) into the stress 
balance equation 

V- < er° >= 0, (A-92) 

gives us three equations for the four unknowns (u° and p ). We now turn 
to continuity to give us the final equation. Continuity will give the final 
equation. Since the fluid stress balance (|A-79j) is linearly forced by the 
external pressure gradient we can define a tensor k relating the external 
pressure gradient to the pore scale flow: 

v° - d t u° = -k • Vp°, (A-93) 

Averaging over the fluid volume yields 

< v° > -<p f d t u° = - < k > -Vp°. (A-94) 

Since the fluid is incompressible averaging the continuity equation 

V • v° + V x , • v 1 = 0, (A-95) 

gives 

= V- < v° > +i J V x , ■ v 1 dV = V- < v° > +i J n • v^S 

= V- < v° > +^ J n • d t u l dS = V-<v°>-^y V - dt^dV, (A-96) 

where we have used l)A-70|l . Taking the divergence of (|A-93|) and using 
(|A-96|1 and (|A-89|) to eliminate v° and u 1 respectively yields 

- V- < k > -Vp° = V • (< v° > -4> f d t vP) 
=< V x , • B >: e(3 t u )- < V x > ■ C > d t p° - d t u° • V<f> f - 0/V • atu^A-g?) 

If the change in solid volume fraction is much smaller than the volume frac- 
tion itself dtu ■ V(j)f ~ 0. Furthermore, if the solid skeleton is incompressible 
then V- < v° >= so that the first two terms on the right hand side of 
equation (|A-97|) are negligible. For a compressible isotropic skeleton (|A-97|) 
yields 

(3d t p° - V- < k > -Vp° = -ad t V ■ u°, (A-98) 

where the (3 =< V x / • C > is the bulk compliance of the solid skeleton and 
a = (pf — < V x ./ • B >a /3 is the effective fluid volume fraction. In general 
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if the solid is treated as compressible, the fluid must also be treated as such 
since their bulk moduli are comparable. Thus, j3 is really a measure of the 
compressibility when the system is jacketed, so that for a mixture of an 
incompressible solid and fluid, (3 = 0. Multiple scale analysis (Auriault h 
Sanchez-Palencia 1977) shows that < V x / ■ B >=< A : e(C) >= 7 so that 
(A27) takes the final form 

< <r° >= 2^e(u°) + AV • u° I - ap°I. (A-99) 

Equations (|A-92|) . QA-98JI and l)A-99|) are the equations of poroelasticity, 
identical in form to the equations written down by Biot (1941). Removing 
the brackets and superscripts we recover equations (JHJ) and (JSJ) from £j21 
Studying poroelasticity from the micro structural point of view allows us 
to see that Biot's (1941) equations correspond to a locally compressible 
solid skeleton and the equations of mixture theory (Barry & Holmes 2001) 
correspond to an incompressible solid skeleton. 

* Appendix B: Equations of motion for a poroelastic sheet 
The equations of motion for a sheet of thickness H and length L are 
found using the same techniques as for a filament. The displacement field 
u = (u(y),v(y),0) is two dimensional, where the y direction is normal to 
the neutral surface and the free surfaces are located at y = ±H/2. We use 
the following dimensionless parameters 

a 2 ^H 2 , 2jxa H 2 , 

+ H 2 , ^ + \)H* 

V = Hy ' axx = ^m? a ™> P = 3(2^ + A)L 2 P ■ (B - 100) 

The dimensionless parameter 5 characterising the ratio of the fluid dtress to 
the solid stress is ^ 

6 = (/x + A)[/3(2^ + A) + a 2 ] (IM01) 
The pressure field is found by solving the 1-dimensional diffusion equation 

d t p-d yy p = -yd xt e, (B-102) 

with the boundary conditions 

dyp + Bip = at y = 1/2, 
-d y p + Bip = at y = -l/2. (B-103) 
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Then 

rt 



P = -Y / Xnsm^y [ e-^-^d^Odt' (B-104) 



where the A n satisfy 



^cos^^ + Bi sin^^ = 0, (B-105) 



and Xn and 7„ are given by 

2(2 + Bi) sin 2(2 + Bi) 2 sin 2 ^ 

Xn = — 1 ^ , 7n= \„, ; ^.^, 2 ■ (B-106) 



An(l-^) A 2 (l 



Equations (|B-105|) and ()B-106|) together with equation (|26|) for the motion 
of a poroelastic plate with time-dependent plane stress. 

*Appendix C: Kirchhoff-Love theory for a bent, twisted filament 
In this section we construct the equilibrium equations for a thin poroe- 
lastic rod whose deformation is not necessarily in the plane. The case of a 
purely elastic filament is treated in Love (1944). The configuration is given 
by the position of the centre line and the orientation of its cross-section 
at every point along it. At every point along the centre line of the rod 
r(X, t) = (X(x,t),Y(x,t), Z(x,t)), where x is the arc-length we consider 
the orthogonal triad dj(x, t), i = 1, 2, 3, where di and d2 lie along the prin- 
cipal axes of the cross-section of the rod and 

d 3 = d x r (C-107) 

is the vector tangent to the centre line. The orientation is determined by 
a body-fixed director frame that allows us to consider finite deformations. 
In this case we use the director basis to follow the evolution of the fluid 
pressure field. The vector of strains k is given by 

K = K (1) d!+K (2) d2 + ftd 3 , (C-108) 

which defines the rotation of the principal axes along the filament. Here 
«W and are the projections of the curvature of the centre-line onto the 
principal axes of the cross-section and Q is the twist strain. 

d x di = Kxd t . (C-109) 

The stress resultant vector F(x,t) and the couple resultant vector M(x,i) 
at any cross-section can be written as 

3 3 

F = ^2F®(x,t)di(x,t), M = ^ (x, t)di (x,t), (C-110) 

1=1 i=l 
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where and are the shear forces and and are the bending 
moments along the principal axes, F& is the tensile force and is the 
twisting moment. 

Since the equation for the diffusion of pressure is linear we consider 
the bending about the principal axes separately. In light of equation (|25|) 
above (d x 9 being the curvature along one of the principal axes) we can write 
equations for the dimensionless couple resultant vector M 

n=l J ° 
n=l J ° 



M {3) = cn, (c-111) 



where C(= ^xp^for a circular rod) is the dimensionless torsional rigidity 

(normalised by the bending contribution to AfW), r a is the dimensionless 
twist strain, and the X n are determined by solving equation (|18[). We note 
that the twisting moment has no poroelastic contribution because it is purely 
a shear deformation, and poroelastic effects arise only from volumetric de- 
formations as seen in equation Finally, the local balance of forces and 
torques give the equilibrium equations 

d x F + F 'external = 0, (C-112) 

d x M + d 3 x F = 0, (C-113) 

where F external is the external body force acting on the cross-section. The 
complete set of equations that determine the poroelastic behaviour of a 
filament are (IC-107ll . (l(J-lllll . (l(J-112fl and (IC-1131) . 
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